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ABSTRACT 

We examine the conditions for the revival of the stalled accretion shock in core-collapse supernovae, 
in the context of the neutrino heating mechanism. We combine one dimensional simulations of the 
shock revival process with a derivation of a quasi-stationary approximation, which is both accurate 
and efficient in predicting the flow. In particular, this approach is used to explore how the evolution of 
the system depends on the shock radius, i?s, and velocity. Vs (in addition to other global properties 
of the system). We do so through a phase space analysis of the shock acceleration, as, in the Rs — Vs 
plane, shown to provide quantitative insights into the initiation of runaway expansion and its nature. 
In the particular case of an initially stationary (Vs = 0, as = 0) profile, the prospects for an explosion 
can be reasonably assessed by the initial signs of the partial derivatives of the shock acceleration, in 
analogy to a linear damped/anti-damped oscillator. If das/dRs < 0 and das/dVs > 0, runaway 
expansion will likely occur after several oscillations, while if das/dRs > 0, runaway expansion will 
commence in a non-oscillatory fashion. These two modes of runaway correspond to low and high mass 
accretion rates, respectively. We also use the quasi-stationary approximation to assess the advection- 
to-heating timescale ratio in the gain region, often used as an explosion proxy. Indeed, this ratio does 
tend to ~ I in conjunction with runaway conditions, but neither this unit value nor the specific choice 
of the gain region as a point of reference appear to be distinct conditions in this regard. 

Subject headings: hydrodynamics - shock waves - instabilities - supernovae: general 


1. INTRODUCTION 

The physical mechanism which drives core collapse 
supernovae remains an outstanding problem after sev¬ 
eral decades of research. While there exists clear evi- 
dence that massive stars do explode (|Smarttl [200911 . a 
viable explosion mechanism has not yet been demon- 
str ated robustly by theoretical m e ans (for recent review s. 
see lBurrow ^ (|20I,'llfl : l,Iankal()20I2[ l:jFo glizzo et aT1(|20I5[B . 

Since first p ropos ed by Iwilsoni (|I985 1 and 
iBethe fc WilsonI (|I985[ 1. ’’delayed neutrino heating” 
has generally been considered the most plausible mech- 
anism fo r drivi ng core collapse supernovae (see however, 
iKushniij (120151) for an alternative view). The process 
envisioned for this mechanism is essentially two-staged. 
First, the iron core of the massive star collapses to a 
proto-neutron-star (PNS) which is stabilized by the 
strong nuclear force, creating an accretion shock that 
plows the incoming material which flows inward. This 
shock stalls, due to heavy energy losses in neutrino 
cooling of the shocked accretion layer (and also through 
dissociation of nuclei as they cross the shock front), 
but is eventually revived when heating of this layer by 
neutrinos emitted from the PNS overcomes the energy 
losses and the inertia of the incoming material. The 
fundamental issue regarding this mechanism is how 
(and if!) the competition between neutrino heating 
and cooling, as well as gravity and ram pressure, can 
revive the accretion shock and drive it into an outgoing 
expansion, eventually disrupting the entire envelope of 
the star. 

After three decades of simulations of neutrino heat¬ 


ing following core collapse, the overall picture is still 
a complicated one. There is a broad consensus that 
self-consistent, one dimensional simulations generally 
fail to explode for all bu t the lowest-mass pro genitors 
(|Liebendorfer et all I200II : iKitaura et al.l 120061) . This 
fact has motivated a shift towards two- and subse¬ 
quently three-dimensional (3D) simulations. In these 
simulations, multi-dimensional effects, such as turbu¬ 
lence, convection, rotation and instabilities, come into 
play, and have been demon s trated t o generally aid a n 
explosion (see, e.g.. iCouc hi (EoTl : lEernandezi (fWl : 
iMclson. Janka fc Marekl ()2ni5[ )). On the other hand, 
state of the art 3D simulations have yet to resolve 
the long standing problem: they tend to predict sub- 
energetic explosions (kinetic energies of a few I0f^_er^_or 
less), or no explosion at all (see the discussion in lBurro-^ 
(|20I3l) and references therein). However, the complexity 
of the neutrino-driven mechanism (diverse in both the 
physical processes and the wide range of distances and 
time scales involved) means that fur ther advanc e s may 
yet change this conclusion (see, e.g., iLentz et'all (|2015l) 
for a more favorable outlook). 

This inherent intricacy has also led to an additional 
class of studies: effective simulations and calculations, 
which use simplified assumptions and physics. Such 
works compromise on quantitative accuracy in order 
to facilitate a qualitative understanding and a more 
straightforward parameter survey for identifying the un¬ 
derlying principles which are necessary to generate an 
explosion. A cornerstone of this line of research has been 
to invoke the neutrino luminosity as a free parameter 
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in the simulations (a ’’neutrino light-bulb”), instead of 
generating it self-consistently in the simulations. Unsur¬ 
prisingly, a fiducial increase of the neutrino luminosity 
does lead to an explosion. These analyses gave rise to 
the concept of a ” critical neutrino luminosity”, which is 
the minimal luminosity that drives a runaway explosion 
when all other parameters of the system are predeter- 
mined. The critical lum inosity was first introduced by 
iBurrows fc Gosh^^ (|1993[) , who considered the problem in 
the stationary approximation, identifying as critical the 
minimum luminosity for which no steady-state solution 
exists. 

Stationary models of similar nature were applied 
also in recent works, in search of the origin of the 
critical nature of the problem, as well as a clearer 
condit io n for an explosio n (see iPejclm^^honigson 
(2012); iKeshet fc Balber^ (120121) : iMurnhy fc Dolence 
( 20151) and references therein). For example, under sim¬ 
plifying assumptions such as a neutrinosphere of black 
body temperature = 4r4 MeV and mass density 
Pn = lO^^pii g cm“^, the cri tical luminosity for whic h 
a solution cannot be found is (|Keshet fc Balber^l2012l) . 


shock to an explosion, and their relation to the critical 
neutrino luminosity. Furthermore, the critical neutrino 
luminosity for an explosion initiated by oscillations ap¬ 
pears to be somewhat lower than that predicted by the 
stationary models, and we aim to uncover the reason for 
this trend. 

The structure of this paper is as follows. In section [2] 
we review the spherically-symmetric physical model used 
in this work. Our simulation code, developed originally 
for this work, is described in section [31 Typical results 
concerning the accretion flow which transitions from os¬ 
cillations to an explosion are shown in section Our 
fundamental observation, considering the conditions for 
a positive shock acceleration, is presented in section |5l 
Here, with the aid of appendices Kl and iBl we develop a 
quasi-stationary approximation which allows us to study 
the shock properties in the phase-space of the shock ra¬ 
dius and velocity. We compare our conclusions to a fre¬ 
quently suggested timescale criterion for an explosion in 
section^ In section|7]we summarize our conclusions and 
discuss some of their implications. The oscillation period 
is estimated in appendix ICl 


Lr ~ 6.9 X 10 


52 
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MiT, 


37^ 


erg s \ (1) 


where Mi = M/{1 Mq s ^) and Mp^i.3 = Mp/1.SMq 
are the normalized mass accretion rate and mass of the 
PNS. 

Treating the neutrino luminosity as a free param¬ 
eter in simulations has also led to the key observa¬ 
tion that the critical luminosity in multi-dimensional 
simulations is lower than in one dimensional simula¬ 
tions under similar condi tions ([Murphv fc Burro'^120081 
iCouch fc O’ conno HMI), thus explicitly highlighting the 
importance of multi-dimensional effects. We note that 
there are conflicting results regarding whether the criti¬ 
cal luminosit y in three dimensions is higher than i n two 
dimensions dBruenn et al.l I2009I iNordhaus et alJ 2010; 


Dolence et al. I 120131 Takiwaki. Kotake fc Suwal 2014 : 
Couch fc O’connoill2014l) . 


In this work we revisit the critical nature of the 
transition from a steady accretion to a runaway explo¬ 
sion, with the aid of effective, one dimensional simula¬ 
tions. While some important multi-dimensional features 
are necessarily discarded when using this approach, it 
still qualitatively describes much of the dynamics which 
dictates the evolution of the accretion shock and the 
shocked material. By nature, one-dimensional simula¬ 
tions are better suited for parameter surveys, being eas¬ 
ier to evaluate as a quantitatively well-defined prob¬ 
lem. We are specifically motivated by the fact that 
in simulations, explosions can occur after the shock 
goes through a series of increasingly strong oscillations, 
rather than directly ac celerating from a standing shock 
to runaway expansion (TO nishi. Kotake fc Yamada 120061 
IMurphv fc BurrowsI 120081 iFernandezl 120121) : naturally, 
this feature cannot be assessed by purely stationary mod¬ 
els (in which the stability of the solution can only crudely 
be investigated (iNagakura et al.l l2013D b Our general 
goal is to examine which quantitative aspects of the flow 
determine the transition from an oscillating accretion 


2. THE SPHERICAL MODEL 


Simulations generally indicate that following core 
bounce, there is a transient phase of order 100 ms, af¬ 
ter which the incoming mass accretion rate, Mq, and 
the neutrino luminosity, L, settle o n roughly fixed values 
during the evolution o f the shock (IBurrows et al.ll2007bl 
iMarek fc Jankal 120091) . Correspondingly, it is common 
practice in simplified models of the explosion process to 
set these two parameters as constants, and study the 
dynamical dependence on their values. Here, we re¬ 
duce the problem to an idealized spherically symmetric 
flow, as has been done in many similar studies (see, e.g., 
iFernand^ ()2012D ). In spherical symmetry, the equations 
of motion, used to calculate the dynamics, are the con¬ 
servation of mass, momentum and energy: 


dp I d , 2 \ n 
717+ 32^1^ P«)=0, 


dt 


and 
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= pq 


( 2 ) 

( 3 ) 

( 4 ) 


where u, p, and p are respectively the fluid radial velocity 
(in the lab, i.e. the PNS, rest frame), mass density, and 
pressure, and 


^tot — 2 ^ 


e — 


GMp 

r 


( 5 ) 


is the specific total energy, and e is the specific internal 
energy. The net neutrino deposition rate (heating minus 
cooling) per unit fluid mass is denoted by q. 

In the gravitational terms (G being Newton’s constant) 
Mp is the PNS mass. In the general case, Mp should 
be replaced with the mass enclosed within a radius r at 
time t, but since in realistic scenarios the mass of the 
PNS dominates over the mass of the accretion layer, we 
neglect this layer’s self-gravity and use a fixed central 
mass. 
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2 . 1 . The Basic Physical Model 

Equations dm are solved for an accretion layer lo¬ 
cated between the PNS and the shock radius, Rs- The 
immediate downstream of the shock is thus used as an 
outer boundary. It is common practice in simplified sim¬ 
ulations of neutrino driven explosions to also specify an 
inner boundary at a radius identified with the PNS sur¬ 
face, Rp. Physically speaking, the radius Rp roughly 
corresponds to the neutrinosphere, above which the neu¬ 
trino luminosity is approximately constant. 

The outer boundary conditions at Rp are determined 
by approximating that the preshocked material as in a 
pressure-less free-fall with a constant mass accretion rate 
Mq = AttR^poUq. The velocity (in the lab frame, i^ the 
PNS frame) and density of infalling material at the shock 
are then 


uq = —a 


(GMpV^^ 
\^) ' 


Po = 


\Mo\ 


47ra(GMp)i/2 


R 


-3/2 


( 6 ) 


where a is a constant of order unity (v^ for perfect free- 
fall). The properties of the material in the immediate 
downstream of the shock are then determined by the 
Rankine-Hugoniot jump conditions. 


PlVi — PqVq , 

Pi + Pivl =P0 + PqvI , 

and 

1 2 A 1 2 Po 

-Vi + Cl -— -Uq -I- Co H- Qd ■ 

2 pi 2 Po 


( 7 ) 

( 8 ) 

( 9 ) 


Here indices 0 and 1 denote quantities upstream and 
downstream of the shock, respectively, and v = u — Vg 
is the fluid velocity relative to the shock, taking into ac¬ 
count the (lab frame) shock velocity, Vg- Hereafter we 
denote pi = pg, ui = ug and Pi = Pg, indicating the 
mass density, fluid (lab frame) velocity, and pressure just 
beneath the shock. 

Finally, qd in equation ([9|) is the energy removed per 
unit mass through the dissociation of infalling ions by 
the shock. The value of qd has several important conse¬ 
quences for the entire profile of the accretion layer. In 
terms of the boundary conditions, it determines the com¬ 
pression ratio across the shock, /3 = pi/po, given by 


7-1-1 


7 - + 2(72 - 1)0 


where 0 = 


Qd 


{uo -Vgf 

( 10 ) 


Here, 7 = P/{pe) -I- 1 is the effective adiabatic in¬ 
dex of the shocked material , which is typ ically radiation- 
dominated and so 7 ^ 1.4 (|Jankal[200lll . In the limit of 
zero dissociation [qd —>■ 0 , 0 —>■ 0 ), the compression factor 
equals that of a simple strong shock, /3—>• ( 7 -|-l)/( 7 — 1 ), 
but if dissociation is significant (compared to the kinetic 
energy of the infalling material in the shock reference 
frame) the compression factor will be larger. 

2 . 2 . Further Simplifications 

As our goal in this work is a qualitative interpretation 
of the conditions for shock revival through neutrino heat¬ 
ing, we apply some further simplifications to our model, 


which allow for a clearer insight into the numerical sim¬ 
ulations. First, we use a simplified equation of state 
for the shocked material, describing radiation, nonrela- 
tivistic nucleons, and relativi stic electrons of zero chemi¬ 
cal potential and degene racy (|Yamasaki fc Yamadall2005l : 
iKeshet fc Balber5l2ni2[l : 


11 ^4 


kppT 


11 


e = 


SkpT 

2 


( 11 ) 


where a, kp, and rrin are the radiation constant, Boltz¬ 
mann constant, and nucleon mass, respectively. While 
this equation of state neglects the finer aspects of the 
composition of the shocked material, it allows for more ef¬ 
ficient simulations while still capturing the main essence 
of the problem, especially the transition from a radiation 
dominated state (7 ~ 4/3) near the shock to a matter 
dominated one (7 ~ 5/3) near the PNS, as the density 
increases by several orders of magnitude. Correspond¬ 
ingly, the dissociation energy, qd, is treated as a free pa¬ 
rameter, which can range between the non-dissociation 
limit qd = 0 and full dissociation of iron ions into free 
nucleons with qd = qpe = 8.5 x 10^® erg g“^. In reality, 
the actual dissociation across the shock is partial, due 
to recombination processes of nucleons into a particles, 
so that some of the dissociation energy is later added 
back to the accreted material. Notice that when dissoci¬ 
ation is included, there is a physical upper limit on the 
radius of the stagnation shock for free-falling material, 
Ps(t = 0 ) < 200Mpp,3{qd/qFe)~^ km. 

We also use a simple recipe for neutrino heating and 
cooling, which is often applied in simplified analytical 
and n umerical models (|MurDhv fc Burro^ 120081 iJankal 
1200111 . The simple formula for the total neutrino heating, 
qn and cooling, qc, are 


qn 


= 1.54 X 102°L52 



and 


4MeV 


2 

ergg“^ s“^ 

( 12 ) 


qc = 1.40 X 10^° 


T 

2MeV 


6 

ergg-^s”^ 


(13) 


where L 52 is the electron-neutrino luminosity in units 
of lO^^erg s“^, is the electron-neutrino temperature 
at the neutrinospehre, and T is the (radius-dependent) 
temperature of the matter and photons, assumed to be in 
equilibrium. The total energy deposition rate in equation 
([3]) is then q = qp — qc. In equation (fT^ it is assumed 
that the electron antineutrino luminosity is equal to the 
electron neutrino luminosity, and that the contribution 
of the other neutrino types to heating can be neglected. 
The total energy deposition rate in equation (jd]) is then 

9 = - qc- 

Finally, since we are only considering matter above the 
neutrinosphere, where the optical depth is small, we as¬ 
sume a fixed neutrino luminosity above radius Rp, so 
that L in Equation m is independent of radius. In 
principle, the heating term can be corrected by a factor 
of e“'^, where r(r) is the optical depth between radius Rp 
and r, but we find that applying this correction has a mi¬ 
nor effect on the bulk properties of the accretion layer, 
which are the focus of our present work. Therefore, we 
neglect this correction. 
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3. THE SIMULATION CODE 

We solve the flow equations (lai) with the aforemen¬ 
tioned boundary conditions and simplifications using 
a one-dimensional Lagrangian code. The code imple¬ 
ments a standard vo n Neumann and Richtmyer stag¬ 
gered mesh method (rron Neumann fc Richtmverl 119501 : 
iRichtmver fc Morton 1 119671 1 for the equations of mo¬ 

tion. The energy equation is solved implicitly with non- 
adiabatic contributions given by equations (dm). 

The outer boundary is followed by continuously adding 
Lagrangian cells above the shock in free-fall velocity ac¬ 
cording to equation ([ 6 ]). The shock dynamics are cal¬ 
culated with a quadratic artificial-viscosity, tr, added to 
the pressure term of the flow, but this requires some ad¬ 
ditional caution. Energy losses by dissociation at the 
shock (when included) can amount to a significant frac¬ 
tion of the internal energy of the post shock material, and 
as a result, numerical disturbances may arise within the 
artificial-viscosity scheme in which the shock is smoothed 
over a few grid cells. We circumvent this problem by ap¬ 
plying the energy loss gradually as the cell passes through 
the shock region. For every cell entering the shock (with 
a > 0.5P), the density and internal energy are saved, 
and the asymptotic post-shock compression factor, /?, is 
calculated according to the density profile of the flow. 
We regulate the amount of energy lost by a cell, qd, so 
that it reaches qd only after the cell is compressed to the 
asymptotic value. Quantitatively, qd is calculated by 

{ 1 if p > Pi; 

if po < p < Pi; (14) 

0 if P < Po , 

where po is the pre-shock density of the element, and 
Pi = /3po is the expected density after the shock. The cell 
then loses up to qd energy, but only as long as the internal 
energy does not drop below its value before entering the 
shock. This recipe guarantees that the internal energy 
in a cell cannot drop to negative values through rapid 
dissociation losses. 

In order to ensure stability in the shock downstream, 
non-adiabatic processes (neutrino cooling and heating) 
are incorporated only after the cell has lost a total en¬ 
ergy of qd- This gradual dissociation loss recipe over a 
typical shock width (a few cells) is a minor correction 
since neutrino heating and cooling are generally much 
weaker near the shock than farther downstream, and we 
find that applying changes the critical luminosity by no 
more than a few percents, while guarantying numerical 
stability. 

At the inner boundary, mass elements that enter the 
PNS (having m < 0 at Rp) interact with a constant pres¬ 
sure Pp in a ghost cell just below Rp, until they drop 
completely below Rp (are absorbed in the PNS) or al¬ 
ternatively attain a positive velocity. The latter occurs 
only when the flow is well into a runaway expansion. 

The simulations were typically calculated with about 
500-1000 cells in the accretion layer, maintaining a de¬ 
creasing cell width toward the PNS: the cell widths 
are adjusted so that each cell is thinner by a factor 
of (1 -I- A)“^ with respect to its upper neighbor, with 
A Ri 10“^ — 10“^. Rezoning is applied when necessary. 
We find that this resolution allows for numerical conver¬ 


gence, both in terms of the flow near the shock and in 
the steep density gradient near the PNS. 

4. FEATURES IN AN EXPLOSIVE FLOW 

We initialize a simulation for a given combination of 
L, Mp, Rp and Mg by determining a stationary pro¬ 
file, i.e, solving the flow equations when all par¬ 

tial time derivatives are set to zero, as is the shock ve¬ 
locity, Vsp = 0) = 0. The outer boundary conditions 
follow by setting a specific shock radius, Rs{t = 0), so 
we may solve the entire profile by integrating the sta¬ 
tionary flow equations from the shock radius inward to 
Rp. The resulting pressure, Pp in the ghost simulation 
cell just below Rp, is then determined, and subsequently 
serves as a time-independent inner boundary condition 
aX Rp. As mentioned above, this pressure is assumed to 
be a property of the (unsimulated) PNS, and is thus kept 
fixed during the simulation. The dynamical evolution is 
then initiated with some small perturbation, either by 
intentionally shifting the shock radius (typically by few 
hundred meters), or simply by numerical noise. We hnd 
that as long as the initial perturbation is small, the evo¬ 
lution that follows does not depend on the specifics. 

We mostly varied L and Mg while keeping the PNS 
mass fixed at Mp = 1.3 Mq, and its radius at Rp = 
40 km. We also varied the dissociation parameter (qd), 
between the non-dissociation limit qd = 0 and full disso¬ 
ciation of iron with qd = qpe- The choice of the initial ra¬ 
dius of the stalled shock Rs{t = 0) warrants some discus¬ 
sion. Without enforcing some additional constraint, this 
radi us has no unique value. In the stationary approxima¬ 

tion (iBu rrows fc (rOsh ^flO^lPeicha fc ThomDSonll2012t 
iKeshet fc Balberdl2012h , the shock radius is uniquely de¬ 

termined, usually by requiring that the optical depth 
between t he PNS and the shock b e equal to 2/3. As 
shown by IKeshet fc Balberd (|2012f) . for such an optical 
depth the shock radius is a slowly varying function of 
L (when other parameters are kept fixed), and is typ¬ 
ically (100 — 250) km. In the dynamic simulations we 
follow here, the accretion layer goes through a range of 
shock radii and optical depths, and so we opt to begin all 
simulations in a specific parameter survey with the same 
initial shock radius, typically Rs{t = 0) = 100 km. In 
T5.4I we vary the initial shock radius, in order to demon¬ 
strate that it too affects the critical neutrino luminosity. 

One of our goals is to identify the process of shock 
revival through the growth of oscillations, u ntil the on¬ 
set of runaway expansion, which - following iFernand^ 
- we define the onset of an oscillatory explosive 
flow. Examples of the time-dependent evolution of the 
shock radius and velocity for Mp = I.SMq, Rp = 40km, 

I Mg I = O.SMq s“^ and different neutrino luminosities are 
presented in Figure [T] Dissociation energy losses at the 
shock were neglected {qd = 0) in these simulations. All 
these simulations are initiated with the accretion shock 
set at Rs{t = 0) = 100 km, with a stationary profile that 
corresponds to this radius. 

As seen in the figure, most simulations show some ini¬ 
tial oscillations of the shock radius and velocity, which 
have typical periods of tens of milliseconds (see appendix 
E] for the reason for this narrow range of periods). The 
distinction between an explosive and a nonexplosive case 
is clearly evident: when the neutrino luminosity is large 
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enough, the oscillations grow in amplitude until eventu¬ 
ally the shock velocity no longer reverts to a negative 
value, but rather continues to increase with a positive 
acceleration, culminating in a runaway expansion. The 
lowest neutrino luminosity which drives such a flow cor¬ 
responds to the critical luminosity - in this case, about 
Lcrit ,52 ~ 8 , while for higher luminosities the onset 
of runaway expansion occurs earlier (after fewer oscil¬ 
lations); for L 52 = 10 , barely one oscillation is com¬ 
pleted before the shock velocity growth becomes expo- 
nential (thi s can b e considered a ” non-oscillatory mode” 
(|Fernanded I2012f) l. The curve for L 52 = 7.5 repre¬ 
sents results typical of L < Lcrit- small oscillations are 
damped and the flow settles back to stationary accretion. 



Figure 1. Simulations of shock dynamics [radius (top panel) 
and velocities (bottom panel)] for various luminosities [L 52 is the 
neutrino luminosity in units of lO^^erg s“^], both bellow (black 
dashed curve) and above (other curves) the critical luminosity. 
Fixed parameters in the simulations are: |Mo| = O.SMq s~^, 
Rs{t = 0 ) = 100 km and = 0 . 



time [s] 


Figure 2. Shock dynamics as in Figure [T] but including dissoci¬ 
ation losses qd = qp^ = 8.5 X lO^^erg g“^. 


This is to be expected, of course, since increased lumi¬ 
nosity is required to compensate for the energy loss rate 
\M\qd fv 1.4 X 10®^ erg s“^. In practice, the critical lu¬ 
minosity is increased by a larger margin than \M\qd, to 
Lcrit ,52 ~ 13.75, since dissociation also changes the com¬ 
pression factor at the shock, and hence the outer bound¬ 
ary conditions. 

5. SHOCK ACCELERATION IN A PHASE-SPACE DIAGRAM 


In this section we develop the quasi-stationary approxi¬ 
mation to derive an approximate expression for the shock 
acceleration, ag, given the shock radius, Rs and veloc¬ 
ity, Vs- This derivation allows us to identify regions in 
the Rs — Vs plane in which both the velocity and ac¬ 
celeration remain invariably positive, hence indicating in 
a phase-space fashion when the conditions are favorable 
for a runaway expansion. 

First, consider the spherical (with respect to the cen¬ 
ter of the PNS) moment of inertia of the accretion layer 
between Rp and Rs, 




Anr'^pdr . 


Its first time derivative is 


dl 

dt 



d'lrr'^-^^dr 

dt 






d-Kr'^^dr + diTRsPsVs ■ 


(15) 


(16) 


The notation [• • • ]p stands for the difference between the 
expression in the immediate downstream of the shock 
and just outside the PNS, where the latter term is null 
in our model since we fix dRp/dt = 0. 

We find that a general feature in the simulated flows 
is that to a very good approximation, the local mass 
accretion rate is uniform in the accretion layer. 


M('r) ~ M = dirr^pu , (17) 


corresponding to slow changes in the density profile. This 
reflects the highly subsonic nature of the downstream 
flow. Note that M is not identical to the incoming mass 
accretion rate, Mq, since it is modified by the shock ve¬ 
locity. Indeed, the mass accretion rate just below the 
shock is then 

Ms = Mo + dnRlVspoid - 1) ■ (18) 

For a radius-independent mass accretion rate, the first 
term on the right hand side of Equation (fT 6 )l is approxi¬ 
mately zero, and so 

L ~ AttRsPsVs . (19) 


In Figure Owe show the results of simulations similar 
to those of Figure [U except that full dissociation losses 
at the shock are included with q^ = qp^. Evidently, 
the inclusion of dissociation losses does not qualitatively 
change the character of the flow; we use this feature in 
when we apply the quasi-stationary model. The in¬ 
clusion of dissociation energy losses does cause a shift in 
the dependence of the flow on the neutrino luminosity, 
requiring higher luminosities to achieve a similar flow. 


Given that dl/dt is now approximated as a function of 
Rs and Vs, the second time derivative of the moment of 
inertia can be used to create an implicit relation for the 
shock acceleration. 





d 



d 

dVs 



( 20 ) 


We emphasize that this is only an approximate equality 
due to the assumption of a uniform accretion rate in the 
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entire accretion layer, resulting in dljdt being a function 
of Rs and Vs alone. 

Equations (I19II20I) relate bulk properties of the accre¬ 
tion layer, dl/dt and d'^I/dt'^, to the quantities at the 
accretion shock. We find that this relationship is re¬ 
produced to a very high accuracy in our simulations. In 
principle, this relation can be used to calculate the shock 
acceleration for a given profile, allowing us to predict the 
oscillatory movement of the shock, and whether or not 
runaway is to be expected. Generally, however, d?I/dt^ 
cannot be derived without a dynamical simulation, for 
which a prediction of as is redundant. Nevertheless, us¬ 
ing some dedicated approximations, we derive in the fol¬ 
lowing a quasi-stationary model, which does allow for the 
prediction oi as, mapping it as a function of the instan¬ 
taneous shock and flow properties. 


5.1. A Quasi-Stationary Approximation 

We modify the simple, stationary approximation by 
including a non-zero shock velocity when setting up the 
conditions at the shock and then solving equations (mign 
as described above. This approximation, which we refer 
to as ’’quasi-stationary”, is applicable when the shock 
velocity is small with respect to the velocities in the ac¬ 
cretion layer, which in turn are smaller than the typical 
sound speed in the flow. Note that the assumption of 
a uniform accretion rate implies that the velocity pro¬ 
file in the accretion layer adjusts quickly to changes in 
the shock radius, while a subsonic flow (also assumed 
in the stationary approximation) ensures that thermo¬ 
dynamic changes in the shock quickly advect throughout 
the shocked material, and influence the inner regions of 
the flow. 

Quantitatively, these conditions can be assessed 
through the typical time scales in the flow, defined by 
the oscillation period, tosc, 


and 


t 


SC 



^adv — 



^env 

M 


( 21 ) 


( 22 ) 


which are the sound crossing time and advection time 
respectively. In the second equality for the advection 
time, is the mass of the envelope, I^ the accretion 

layer. The second equality for tadv is exact only if the 
mass accretion rate is indeed uniform in the entire layer. 
To be precise, tosc is relevant when the flow goes through 
several oscillations before growing to runaway expansion 
(or damping out). Once the flow evolves to runaway 
expansion — or alternatively, if the expansion is initially 
non-oscillatory — the appropriate measure of the shock 
evolution should be its dynamical time, ts = Rs/Vs- 
In the parameter range of interest the post shock flow 
is very subsonic, with tsc being a few milliseconds and 
the shortest of the three timescales. The competition be¬ 
tween the oscillation and advection time scales is more 
complicated, since the latter is very sensitive to the com¬ 
pression factor at the shock, and hence to dissociation 
losses. When neglecting dissociation losses we find that 
for luminosities close to or exceeding Lgrit, the advec¬ 
tion time scale is tadv ~ 10 ms, while the oscillation pe¬ 


riod is about 50 — 60 ms. This typical value for the os¬ 
cillation period can be quantitatively assessed with our 
quasi-stationary model; see Appendix [C] 

The hierarchy tgc < tadv < tosc implies that the quasi- 
stationary approximation should be applicable. We 
demonstrate this explicitly in Figure [21 which compares 
four snapshots of the velocity, specific energy, and mass 
density profiles in the accretion layer, for the case |Mo| = 
O. 8 M 0 s“^, L 52 = 8 (the solid blue curve in Figure [J), as 
found in the simulation and in the quasi-stationary ap¬ 
proximation (using the instantaneous shock radius and 
velocity from the simulation). We also show the sta¬ 
tionary profiles (found for the instantaneous shock ra¬ 
dius but setting Vs = 0) at every snapshot. Clearly, the 
quasi-stationary approximation offers a significant im¬ 
provement over the stationary case, offering an almost 
exact fit not only in specific energy and density, but also 
in the velocity profile of the flow. Correspondingly, the 
quasi-stationary approximation we use below to analyze 
the shock acceleration is suitable for the regime of small 
dissociation losses - or gradual dissociation in the context 
of the microphysics of the shocked material. 

The quasi-stationary approximation is not quite as suc¬ 
cessful in the opposite limit of full dissociation into free 
nucleons across the shock. For qd = qpe, the compres¬ 
sion factor at the shock is higher, generally between 10 
and 20, which results in a low fluid velocity; for lumi¬ 
nosities close to critical we find that tadv can reach 20 
ms. On the other hand, toscS is relatively insensitive to 
the level of dissociation, and so the quantitative devia¬ 
tions of the actual shock acceleration from that predicted 
by the quasi-stationary model become larger. Nonethe¬ 
less, we find that the qualitative behaviour of the flow 
and runaway expansion is similar even in the limit of 
full dissociation (since the oscillation period remains the 
largest time scale in the problem), and so we do assess 
that the quasi-stationary approximation is generally a 
useful starting point in the qualitative analysis of the 
transition to runaway expansion. 

The advantage of the quasi-stationary approximation 
is that it allows us to estimate the second time derivative 
of the moment of inertia, [d^I/dt^)Qs, which has units 
of energy. The full derivation is given in Appendix 0 
here we quote the final result, 

1 d^/ 

+ Wb , (23) 

where the three components on the right hand side are 
an effective energy, Eq 5 , a work term associated with 
the PNS, WpNS, and the energy advected across both 
boundaries, Wb- These are defined as follows. 

The effective energy, composed of kinetic, gravita¬ 
tional, and internal contributions, can be written in the 
form 

v^Qs-k + h + jj + Bs-k + h + u, (24) 

where the last expression can be estimated from the 
quasi-steady model. Here, we defined 
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t=0.5 [s] t=0.6 [s] t=0.7 [s] t=0.8 [s] 



radius [km] 


Figure 3. Spherical symmetric simulation compared to quasi-static (QS) approximation and steady-state (SS) approximation, for \Mq\ = 
O.SMq s~^, L 52 = 8 , Rs{t = 0) = 100 km and = 0. Shown are the velocity (top), specific energy density (middle), and mass density 
(bottom), at four different times (see labels). The smoothing of t he d ensity profile at the simulated shock is a numerical effect, due to the 
use of artificial viscosity and the recipe for dissociation (equation [Till. 


and 




dm , 




6 s(m) dm , 


a moving shock, are isolated in Bg and subsequently ne- 
' ' glected; see appendix [A] 

Next, there is work done by the PNS, arising from 
static and ram pressures at i?p. In general, static pres¬ 
sure dominates, so that WpNS can be approximated by 
(28) (see Appendix IBl): 


where 7 is the local adiabatic index, Q is the rate of 
change of the total internal energy due to heating and 
cooling in a layer extending between r and the shock, 


WpNS — 47ri?pPp . (31) 

Finally, the energy advected across the boundaries is 


rRs pRs 

Q{r) = q dm = / Anpr^q dr , 
Jr Jr 


(29) 


and 6 s is the Bernoulli function (specihc energy), 


Wb ^ a(GMp)i/2|Mo|i?y^ 
|Mo| d 


7a{GMpy/^ dt 


(/3-1) 



(32) 



7 _ ^ 2 

6 = -u 


P 

P 


GAA 


(30) 


evaluated for every mass element according to its value 
when crossing the shock. The non-adiabatic (heating 
and cooling) contribution to the energy is contained in 
the double integral U, which tracks the total gained non- 
adiabatic energy in the accretion layer. 

The non-standard coefficients of the energy integrands 
in Eqs. arise from writing Egs in a form that 

is susceptible to the quasi-steady approximation. In this 
approximation, terms associated with the instantaneous 
advection that is implied by assuming a steady state with 


where the second term is the correction that arises due 
to the finite shock velocity. 

Substituting our result for {d^I/dt^Qs in Equation 
((2ni) . we arrive at a closed form expression for the shock 
acceleration in the quasi-stationary approximation: 


as 


Eqs + WpNS - \MoRsuo\ - Vg-^ 

xC 


(33) 


where 

{P-l){Rl-Ry)\My 

2|mo| 


(34) 






























































and 


All the variables in equation (1551) are determined by 
the set of external parameters of the flow (including 
Pp, which we fix at the beginning of the simulation), 
and the dynamical variables Rs and Vs (recall that 
P = P{Rs,Vs,Mp,Mo)). 

5.2. The asiRs^ys) Phase Space 

With the aid of the quasi-stationary approximation, 
we calculate the partition of the Rs — Vs plane into re¬ 
gions of positive and negative acceleration of the shock. 
Given a set of external parameters, Mq, L, Mp and Rp, 
and the conditions of the initially stationary shock, de¬ 
termined by Rs{t = 0) and Vs(t = 0) = 0 (which also 
dictate the boundary conditions at the PNS), every point 
in the phase-space is calculated with the relevant quasi- 
stationary profile, yielding an estimate of the shock ac¬ 
celeration. A set of examples with \Mo\ = O.SMq s“^, 
Mp = I.SMq, Rp = 40 km and Rs(t = 0) = 100 km 
are shown in Figure S) The figure shows the predicted 
values of as in the Rs — Vs plane for three neutrino lu¬ 
minosities which are sufficient to drive a runaway expan¬ 
sion (see Figure [T|). The actual evolutionary path of the 
shock radius and velocity, as found in each simulation, is 
superimposed on each of the plots. 

The most prominent feature in Figure [4] is the distinct 
separation of the Rs — Vs plane into regions of positive 
and negative acceleration. At small and large shock radii 
the acceleration is positive, separated by a trough of neg¬ 
ative acceleration. In the region of small shock radii, the 
acceleration is dominated by the effect of the PNS sur¬ 
face pressure, which is why in this region the dependence 
on the shock velocity is relatively weak, while the spatial 
gradient of the acceleration is quite steep. In contrast, 
for larger shock radii, a positive acceleration depends on 
a more complicated combination of the total energy in 
the flow, and on the boundary conditions at the shock 
itself, both of which depend on the shock velocity as well. 

A key result is the fact that the actual evolutionary 
paths of the shocks in the Rs — Vs plane follow the pre¬ 
dicted acceleration of the quasi-stationary model very 
closely, oscillating between negative and positive veloci¬ 
ties in accordance with the sign of the acceleration. We 
find that the approximation for the shock acceleration is 
qualitatively robust, but tends to overestimate its magni¬ 
tude when compared to the full simulations (the term Bs 
neglected above does tend to reduce |as|). Nonetheless, 
the contours of zero acceleration usually correspond very 
accurately to the points along the path where the shock 
acceleration in the simulations changes sign, especially 
for oscillating shocks. This result lends strong support 
to the validity of the quasi-stationary approach and the 
phase-space analysis. 

The transition to a runaway expansion occurs when the 
oscillations grow sufficiently large such that a velocity of 
order 1000 km s“^ pushes the shock through the negative 
acceleration trough into the region of positive accelera¬ 
tion on the right hand side of the phase space. Typically, 
the path passes close to the saddle point of as{Rs,Vs) 
and arrives in the right-handed region of positive accel¬ 
eration with a positive velocity, so the shock continues 


to expand. Hence, when the shock arrives in this region, 
the oscillations (if present) cease, and exponential expan¬ 
sion ensues. We hereon refer to the right curve (larger 
Rs) of zero acceleration as the ’’critical” as = 0 curve. 
In contrast, the left 05 = 0 curve (smaller Rs) corre¬ 
sponds to the change of sign of the shock acceleration 
in oscillatory motion, and does not necessarily signify a 
transition to runaway expansion. We refer to this curve 
as the ’’oscillatory” as = 0 curve. 

It is important to note that we do not find a case where 
the shock acceleration becomes negative at larger radii, 
to the right of the critical curve. The reason for this 
trend is that as the shock radius increases further, the 
competition between heating and other energetic effects 
maintains a positive acceleration (Equation 1331) . Specif¬ 
ically, while O and {—\MoRsuo\) become more negative 
with increasing shock radius, the effective heating term, 
U, becomes larger (more positive), and generally dom¬ 
inates the change in as (we note that WpNS remains 
constant while K and the last term in Equation (1551) are 
essentially negligible). Hence, once the evolving shock 
crosses the critical as from left to right, runaway expan¬ 
sion is inevitable. 

The magnitude of the neutrino luminosity has a sub¬ 
stantial effect on the acceleration phase-space, as is read¬ 
ily seen when comparing the three plots in Figure |4l 
Clearly, as we raise L the negative acceleration trough 
becomes narrower and shallower, and so the shock re¬ 
quires fewer oscillations to reach a point in the Rs — Vs 
plane which will carry it over to the positive accelera¬ 
tion region. To further elucidate this point, we show in 
Figure [5] the phase space and evolutionary paths for a 
subcritical luminosity of L 52 = 6 (see Figure [1]), for a 
transitional luminosity L 52 = 9.7, and for a high lumi¬ 
nosity of L 52 = 11. In the subcritical case, small per - 
turbations are damped due to local stability (see 115.31) . 
eventually converging onto a stationary solution, whereas 
larger perturbations are bounded due the wide and deep 
negative acceleration area in the phase-space. As the 
luminosity increases, criticality is reached as the damp¬ 
ing turns into anti-damping, and the shock continues 
to expand through a series of oscillations. Increasing 
the luminosity further narrows the negative accelera¬ 
tion trough until eventually the two 05 = 0 curves in¬ 
tersect. At even higher luminosities the curves detach 
again, and this detachment reflects a qualitative change 
in the explosive evolution. Now the initial profile lies 
on the critical as = 0 curve rather than on the oscil¬ 
latory one (recall that we dictate initial profiles with 
Vs{Rs = 100 km) = 0 and as{Rs = 100 km) = 0 ). 
Such an initial profile will necessarily end in runaway ex¬ 
pansion if the initial perturbation drives it into larger 
shock radii, but we find that it may also reach runaway 
expansion after one oscillation, if the perturbation is in 
the opposite, radially negative, direction. We conclude 
that the intersection of the 05 = 0 curves corresponds 
to a transition between oscillatory and direct runaway 
expansion. 

5.3. Paths to Runaway Expansion 

The phase space interpretation leads us to the follow¬ 
ing conclusion: even when all external parameters - {Mq, 
Mp, i?p, Pp and L} - are set, both the initial shock 
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Figure 4. Propagation of the shock in the Rs — Vs phase-space for oscillatory (L 52 = 8,8.5) and marginally non-oscillatory (L 52 = 9) 
expansion (note that we define shocks that runaway after a single cycle as non-oscillatory). The quasi-stationary acceleration, calculated 
in equation EJ, is plotted as contour lines in units of 10® km s Flow parameters are: |Mo| = O. 8 M 0 S Rs{i = 0) = 100 km, and 
qd = 0. Arrows indicate the direction in which the simulated shock evolves, and are marked every 25 ms. 
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Figure 5. Same as Figure!^ but for: top panel - L 52 = 6 (damped oscillations); middle panel - L 52 = 9.7 (near the point where the as = 0 
curves intersect); and lower panel: L 52 = 11 - an initial profile on the critical curve with exponential runaway, either directly (dashed) or 
after one oscillation (solid). 
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radius and velocity determine whether or not the flow 
will become unstable and reach runaway expansion. We 
demonstrate this by means of an example in Figure [ 6 l in 
which we repeat the simulations in the case of L 52 = 7.5, 
including the same pressure at the PNS, but start with 
static Vs{t = 0) = 0 shocks at different shock radii. In 
order to show the evolutionary paths corresponding to 
these different radii superimposed on a single as(i? 5 , Vg) 
phase space map, we do not adjust the boundary PNS 
pressure at according to each initial shock radius {Pp is 
set for Rs{t = 0) = 100 km). Consequently, most pro¬ 
files have a non-zero initial shock acceleration, and evolve 
to cover a wide range of shock radii and velocities. The 
plot shows that there exist small and large initial radii 
for which the flow results in a runaway expansion, while 
in an intermediate regime of shock radii it does not. Con¬ 
versely, we can consider any combination of {Rs, Vs} as 
an initial profile for the simulation. There is a finite re¬ 
gion in the Rs — Vs plane (its border lying between the 
paths labeled Rs = 130 km and Rs = 140 km) for which 
oscillations are eventually damped, and the flow settles 
on the stationary solution. The physical reason for this 
limited range of {Rs, Vs} for which runaway expansion 
is denied (in this particular setup) is that at sufficiently 
small radii the pressure close to the PNS creates a large 
positive shock acceleration, which can drive the shock 
all the way beyond the critical as = 0 curve. Hence, 
runaway expansion will occur either if the flow is initi¬ 
ated in this region, or if is initiated in combinations of 
{i?s, Vs} which eventually evolve so that the shock pen¬ 
etrates deeply into the left handed side region of positive 
acceleration. Of course, the critical as = 0 curve serves 
as an additional limit: a profile which is initiated in that 
region will accelerate exponentially if the initial shock 
velocity is positive (or even only slightly negative). 

In the general case of arbitrary initial {i?S)Vs} the 
outcome of the evolution can be found by calculating the 
entire evolutionary sequence (directly or by following it 
after mapping the appropriate phase space). Moreover, 
when the problem is restricted to initially stationary pro¬ 
files {Vs(t = 0) = 0, asit = 0) = 0), the situation is 
greatly clarified and we find that the outcome of the evo¬ 
lution can be assessed by the initial signs of the partial 
derivatives das/dRs and das/dVs- This can be seen by 
noting that here, perturbations are governed by 

Rs = as = - Rs,o) + -^^^s , (36) 

oRs oVs 

which is simply the damped (or anti-damped) linear os¬ 
cillator equation. This suggests the following behavior. 

• An initially stationary profile with das/dRs > 0 
necessarily corresponds to a point on the critical 
as = 0 curve. As is the case in the lower panel of 
Figure [5l we find that such a profile will invariably 
evolve to a runaway expansion in a non-oscillatory 
fashion. In essence, this class of initially stationary 
profiles is inherently unstable. 

• An initially stationary profile with das/dRs < 0 
necessarily corresponds to a {Rs, Vs = 0 } point on 
the oscillatory as = 0 curve. Whether an initial 
perturbation diverges or damps out can be gauged 
locally by the sign of das/dVs- If das/dVs < 0 , 


we expect a stable scenario: oscillations will tend to 
damp out, and the flow will settle onto the station¬ 
ary solution. If this derivative is positive, runaway 
expansion is expected, but its nature is more sub¬ 
tle. In a purely local perturbation theory, we would 
expect that a further distinction will arise from the 
sign of A = {das/dVsY + '^das/dRs- If A < 0, 
the shock should evolve through a series of increas¬ 
ing (” anti-damped”) oscillations, finally resulting 
in runaway expansion. On the other hand, A > 0, 
along with das/dVs > 0 leads to non-oscillatory, 
runaway expansion (even though das/dRs < 0 ). 

We test this insight in Figure [T] In the figure we show 
the results of a survey of simulations, presenting the crit¬ 
ical luminosities as a function of mass accretion rate for 
Mp = 1.3 Mq and i?p = 40 km. Each simulation is 
initiated with a small perturbation around a stationary 
profile at Rs = 100 km. The survey allows to distinguish 
between stable configurations and runaway expansion, 
and also to distinguish whether runaway occurs through 
an oscillatory or a non-oscillatory path. We define run¬ 
away expansion as non-oscillatory when there is no more 
than one cycle where the shock radius drops below its 
initial value; for example, the L 52 = 9 case in Figure[T]is 
an example of a marginally non-oscillatory runaway ex¬ 
pansion. The figure shows the critical luminosities for os¬ 
cillatory and non-oscillatory runaway expansion, Lc,osci 
and Lc,dir, respectively. These are compared to curves 
which show the loci of das/dRs = 0, das/dVs = 0, and 
A = 0 as calculated with the quasi-stationary approxi¬ 
mation (i.e., finding all combinations of L and \Mq\ for 
which an initial stationary profile with Rs{t = 0) = 100 
km would yield these derivatives). 

As seen in the figure, there is a good qualitative consis¬ 
tency between the simulated critical luminosities and the 
theoretical curves, evaluated by the shock acceleration in 
the quasi-stationary approximation. The das/dRs = 0 
curve is especially indicative of the threshold for runaway 
expansion through the non-oscillatory path at higher 
mass accretion rates, for which das/dVs < 0 for all 
luminosities. Interestingly, the das/dVs = 0 curve 
strongly depends on the accretion rate, requiring exceed¬ 
ingly higher L for increasing \Mq\, eventually diverging 
at \Mo\ — Me = 0.91 Mq s“^. On the other hand, the 
das/dRs = 0 curve has a rather weak dependence on 
the mass accretion rate, so for |Mo| > Me only one class 
of initial stationary profiles exists: this simply reflects 
the fact that for the higher accretion rates, initial pro¬ 
files with Rs{t = 0) = 100 km that explode lie on the 
critical as = 0 curve, rather than the oscillatory one. 
As is indeed confirmed in the simulations, runaway ex¬ 
pansion in this regime can occur only through a non- 
oscillatory evolution. We note that at very high mass 
accretion rates, the as phase space analysis becomes in¬ 
creasingly inaccurate due to a large contribution from the 
Bernoulli function, which cannot be neglected in highly 
non-oscillatory expansions, when the shock dynamical 
time becomes comparable to the advection time (see Ap¬ 
pendix!^. 

The parameter region where das/dRs < 0 displays 
a somewhat more complicated behavior, although also 
qualitatively consistent with the quasi-stationary ap- 
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Rs [k.m] 

Figure 6 . Evolution of initially static (Vs = 0) shocks with different initial radii in the phase-space diagram. In all simulations ^52 = 7.5, 
|Mol = O. 8 M 0 s~^ and = 0. The PNS pressure (and the overall phase space) is calculated for a stable accretion shock at Rgfi = 100 
km. The Flow is unstable for initial radii outside of 84km < Rs{t = 0) < 150 km 


proximation. First, when das/dVs < 0 (low luminosi¬ 
ties), runaway expansion is inherently denied: initial per¬ 
turbations are damped, and the system relaxes to the 
stationary profile, as expected. The das/dVs = 0 curve, 
above which oscillations are expected to anti-damp and 
lead to runaway expansion, indeed lies very close to the 
critical luminosity curve Tc.osc, and hence provides a very 
good estimate for the critical condition for this type of 
evolution. The two curves do not quite coincide, and for 
most Mq there is a small range of luminosities for which 
das/dVs > 0 but runaway expansion does not yet occur. 
We find that in this range oscillations commence, but 
they to tend to settle on a constant or semi-constant (in¬ 
creasing very slowly) amplitude. An example of such an 
evolution is presented in Figure[ 8 l which shows the time- 
dependent radius of the shock for \Mo\ = 0.4Mq s“^ 
and different neutrino luminosities. There is a general 
resemblance to Figure [H but the case of L 52 = 2.5 is 
unique: initially the oscillations grow, but then settle on 
an approximately constant amplitude. Presumably, the 
existence of such an evolutionary path is due to the non¬ 
local nature of the as{Vs) relation (in other words, the 
oscillator in Equation (1361) reaches amplitudes where it 
becomes non-linear). 

Figure [7] also shows the limiting luminosities for which 
there is a transition from oscillatory to non-oscillatory 
runaway expansion in the das/dVs > 0 region. Accord¬ 
ing to linear theory, this curve should coincide with the 
A = 0 curve (which is very close to the das/dRs = 0 
curve, and so is of little importance for the phase space 
analysis). However, in this particular aspect we are 
strongly constrained by numerical resolution, which pre¬ 
vents us from initiating the simulation with truly in¬ 
finitesimal perturbations. In general, it is impractical 
to initiate the simulations with perturbations that are 


smaller in amplitude than about two hundred meters, 
and the eventual dynamics of the shock depend on the 
changes of das/dVs and das/dRs over this scale. In the 
simulations, a perturbation of this magnitude can grow in 
radius by about one kilometer during a single oscillation, 
which is significant in terms of the gradients of the ac¬ 
celeration. We show this by plotting in Figure [7] a curve 
which corresponds to a growth of the initial perturbation 
to one kilometer in a single oscillation Q. The quantita¬ 
tive fit to the critical non-oscillatory luminosity suggests 
that finite resolution leads us to identify non-oscillatory 
runaway when linear analysis breaks down over a 1 km 
scale - obviously below the A = 0 curve. In any case, 
this transitional luminosity is less important than the 
critical luminosities, which separate exploding and non¬ 
exploding profiles. 

5.4. Dependence on the Initial Shock Radius 

So far we have set the initial shock radius arbitrarily at 
Rs{t = 0) = 100 km. Our phase space analysis suggests 
that there should be a non-trivial dependence of the evo¬ 
lution, and therefore of the critical luminosity on this ra¬ 
dius, even when the set of parameters Mq, Mp, and Rp is 
fixed. In Figure[5]we confirm this expectation for the case 
of initial stationary profiles, where we compare the criti¬ 
cal luminosity as a function of the external accretion rate, 
Acrit(Afo), when Mp = 1.3 Mq, Rp = 40 km and the 
initial profiles were calculated for Rs{t = 0) = 100 km 
and Rs{t = 0) = 120 km. We also compare the re¬ 
sults with the critical luminosity found when the initial 
conditions are determined by requiring a r = 2/3 neu- 

^ In a linear oscillator, the amplitude of a perturbation will grow 
by a factor of exp[2'K{das/dVs)/\/\A\\, SO a one kilometer ampli¬ 
tude will be reached in a single oscillation from the initial pertur¬ 
bation when the growth factor is four to five. 
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Figure 7. The critical luminosity as a function of the external mass accretion rate. Single points summarize a survey of simulations and 
show the minimal neutrino luminosities for stable oscillations (dark green squares), oscillatory runaway (green dots), and non-oscillatory 
runaway (red dots). The black line represents the combined critical luminosity. These data are based on our highest resolution, and are 
accurate to within 5%; error bars, which are too small to see at low |Mo|i represent the Icr distribution of the numerical survey. Also 
shown are curves calculated with the quasi-stationary approximation: das/dVs = 0 (orange dashed line; das/dVs > 0 above this curve), 
das/dRs = 0 (blue dashed line; das/dRs > 0 above this curve), A = 0 (gray dashed; A > 0 above this curve), and when perturbations 
are expected to grow to ~ 1km in a single oscillation period (gray dash-dot). See text for details. 
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Figure 8. Same as Figure [T| but for |Mo| = O.4M0 s“^. Note 
that increasing the luminosity changes the evolution from damped 
oscillations to stable oscillations, and only then to oscillatory run¬ 
away, and finally to non-oscillatory runaway. 


trino optical depth through the accretion layer, so that 
Rs{t = 0) is constrained for each simulation separately; 
for all but the lowest \Mo\, this condition forces a larger 
initial shock radius, with Rs{t = 0) = 150—160 km. It is 
evident from the figure that a dependence on Rs{t = 0) 
does indeed exist (weak at small M, becoming more pro¬ 
nounced for larger accretion rates). Our results are con¬ 


sistent with the conclusions of iFernandezi (|2012ll , and are 
distinct from the stationary model, which generates a sin¬ 
gle physical solution when constrained by an additional 
condition. 

Also shown in figure [9] is a similar comparison of sets of 
critical luminosities found in simulations which included 
full dissociation across the shock, qd = 8.5 x 10^® erg g“^. 
The qualitative dependence upon = 0) is retained 
when full dissociation is included, while all the critical lu¬ 
minosities are shifted upwards, (as expected). Notably, 
in this case the shock radius corresponding to r = 2/3 
generally lies between 100 and 120 km, and so the criti¬ 
cal luminosity curve for a constant optical depth is now 
second from the bottom, instead of first . Once again, 
the general trend of the results suggests that our princi¬ 
pal conclusions can be carried over to the realistic case 
when dissociation is accounted for across the shock and 
later d uring accretion. 

The iBurrows fc Gosh"^ (|1993f) critical luminosity, 
Lc,bg{M), is determined by the absence of a stationary 
solution. In cases where a stationary t = solution 
exists, but its instability leads to runaway expansion, the 
critical luminosity in simulations, Lc{M), must be lower 
than Lc^bg{M). A more subtle issue here is that when 
using the r = 2/3 constraint, L^^bg is generally only a 
mild overestimate of Lc even though the issue of stabil- 
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Figure 9. The critical luminosity as a function of the mass ac¬ 
cretion rate for different initial values of the shock radius. Blue 
and red lines correspond to a fixed, Rsi^' = 0) of 100 and 120 
km, respectively, and the green line for an initial shock radius that 
corresponds to an initial optical depth of r = 2/3 in the accretion 
layer. Dashed curves: simulations with no dissociation across the 
shock, Qd = 0; solid curves: simulations with a full dissociation 
across the shock, = QFe- 

ity is ignored (|Fernandedl2012f) . We qualitatively relate 
this to the phase space analysis by noticing that in sta¬ 
tionary models, as the neutrino luminosity is increased 
towards Fc.bg, the initial shock radius changes signif¬ 
icantly - see Figure [TUI As a result, the phase space 
structure of as{Rs, Vs) changes in a substantial manner 
at luminosities close to (but lower than) Lc^bg, trans¬ 
forming the system to a non-stable one. Thus, the criti¬ 
cal luminosities corresponding to the disappearance of a 
solution and to its instability lie in close proximity. 



Figure 10. The initial radius of stationary profiles as a function 
of the neutrino luminosity when constrained by a total neutrino 
optical depth of the accretion layer of t = 2/3. The curves corre¬ 
spond to models with no dissociation across the shock (solid blue 
curve) and full = qpe) dissociation (dashed red curve). The 
curves are ter minated at a luminosity f or which no solution can 
be found (the [Burrows fc GosTl^ 1199311 limit). In these models 
Rp = 40 km and \Mq\ = O.SMq s“^. 


6. A TIMESCALE CONDITION FOR RUNAWAY 
EXPANSION? 

Our phase space interpretation, combined with the 
quasi-stationary approximation, enabled us to identify 
the initial signs of the partial derivatives of the shock ac¬ 
celeration as a criterion for the prospects of an explosion 
when starting from a stationary profile. We emphasize 
that even though they are expressed by quantities at the 


shock, these are global criteria, in the sense that they 
depend on the properties of the entire accretion layer. 
In this section we use the phase space analysis and the 
quasi-stationary approximation to evaluate a commonly 
discussed timescale criterion for the onset of an explo¬ 
sion, comparing heating and advection in the so - called 
’’gain region” of the flow (iMurphv fc BurrowsI l2008t 
Peicha fc ThompsonI l2012t iFernanded l2012t lOtt et ahl 

20m ^ 

The specific neutrino heating scales as the inverse of 
the radial position of the mass element in the accretion 
layer, while the specific neutrino cooling scales as the 
sixth power of the element’s temperature (recall Equa¬ 
tions (11211131) 1. As a result, at smaller radii where the 
temperatures are higher, cooling dominates over heat¬ 
ing, q{r) < 0, whereas closer to the shock heating dom¬ 
inates and q{r) > 0. Such a structure is also found 
in supernovae simulations, including studies which use 
more accurate modeling than ou r simplified equations 
(iFernandezI 120121 : lOtt et al.ll2013D . Correspondingly, it 
is common practice to identify the radial position in the 
flow where q = 0 as the ’’gain radius”, Rq, and refer to 
the region extending between this radius and the shock 
as the ” gain region”. 

The argument concerning the timescales is that if in 
the gain region the heating time. 


^heat,G 



(37) 


is smaller than the corresponding advection time in this 
region. 


iadv,G — 


rRs 

Jrg 


dr 


Mg 
\M\ ’ 


(38) 


then the flow should evolve towards runaway expansion. 
In these equations, Eq and Mq are the total internal 
energy of matter and the total mass in the gain re¬ 
gion, respectively, and the second equality for tadv is ex¬ 
act if the mass accretion rate is uniform. The criterion 
theat,G < tadv.G reflects the expectation that if the ma¬ 
terial can heat significantly before it passes through the 
gain region, steady state accretion could be disturbed. 
In the following, we demonstrate that the phase space 
structure of as can be approximately related to this ra¬ 
tio of timescales. 

First, note that in Equation (l33l) the prefactor C is 
invariably positive, and x is positive for positive shock 
velocities {Vs{P — 1) tends to increase with shock veloc¬ 
ity). Hence, the sign of as is determined by the expres¬ 
sion inside the brackets in this equation. This expression 
includes the effective energy Eqs Ki K + U + which 
characterizes the entire accretion layer, and three addi¬ 
tional terms which include contributions from the PNS 
and the shock. In Appendix iBlwe show that in the region 
of interest, the PNS term Wpns — 4:TTRpPp dominates 
over the advection terms, and so we can generally claim 
that Eq 5 > 0 along with Vg > 0 provide a sufficient con¬ 
dition for runaway expansion. We show this explicitly in 
Eigure (TT] in which we plot again the shock accelera¬ 
tion in the Rs — Vs plane, and superimpose upon it the 
Eqs > 0 curve. As expected, the Eqs > 0 curve lies in 
the runaway expansion region to the right of the critical 
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as = 0 curve, clarifying that a positive effective energy 
is indeed a sufficient - but not necessary - condition. 

The total energy Egg is the result of an integration 
over the entire accretion layer. Conversely, we can define 
a radius-dependent effective energy, E(r), corresponding 
to the integrals (I25II29I) taken from r to Rs- In general, 
we expect E(r) > E(i?p) = Egs, since matter close to 
the PNS experiences the strongest cooling and gravita¬ 
tional potential. As a result, there can exist combina¬ 
tions of Rs and V 5 for which the profile has E(r) > 0 
for some range of radii Rp < r < Rs, while Eq 5 < 0. 
Numerically, we find that for most initial configurations 
that eventually lead to runaway expansion, the function 
E(r) has a maximum at some inner radius in the ac¬ 
cretion layer, which we denote as Rp- This radius is 
always close to, but slightly lower than Rq- just below 
the gain radius the internal energy still increases inward 
due to compression [dU/dr < 0 ), and this effect offsets 
the net cooling by neutrinos. We demonstrate the lo¬ 
cations of Re and Rq in Figure [121 which shows E(r) 
in the accretion layer calculated for various initial shock 
radii and a given set of external parameters (T 52 = 8 , 
\Mo\ = 0.8 Mq Mp = 1.3 Mq, and qd = 0). 

Repeating the process of identifying and calculating 
Re and Rg as a function of Rs and Vs for a given set of 
external conditions, we find the loci of all the points in 
the Rs — Vs plane for which E(i?G) = 0 and E(i?s;) = 0. 
The two corresponding curves are also shown in Figure 
dB and can be seen to lie very close to each other, and 
are both to the left of the critical as = 0 curve. Our 
conclusion is that E{Re) > 0 and E{Rg) > 0 are both 
thresholds that must be crossed when the flow evolves to 
runaway expansion, but they are not sufficient conditions 
conditions to ensure that runaway will occur. Indeed, for 
a wide range of parameters, crossing the E{Rg) = 0 curve 
does not ensure runaway expansion. 

Finally, the quasi-stationary approximation can be 
used to relate the condition tadv,G/theat,G > 1 to the 
condit ion E{Rg ) > 0. With some manipulation of Equa¬ 
tions (IA14llA2f]) . we derive that 

E^Rg) ^ Kg + Ug + i2^Ea - ^) Eg (39) 

with Kg, Eg and Ug being the total kinetic energy, total 
internal energy and energy-gain between Rg and Rs, 
respectiv ely (t he latter is the gain region equivalent to 
equation IA15I) . and ^Eg is the mean adiabatic index of 
the material in the gain region, weighted by the internal 
energy: 

1 

= — / qedm . (40) 

Eg Jrg 

Since Eg is positive by construction, the condition for 
E(i?G) > 0 becomes: 


Eg 


1 

Eg 



|^Q(r)dm > 3 - 2^Ec 


(41) 


where 77 is a numerical factor reflecting the gradient of 
Q{r) (77 = 1/2 for a linear function). We find that for the 
typical configurations in Figure IH 77 ~ 0.6 and ■^Eg ~ 1-4. 

Finally, since the accretion layer typically involves 
small velocities, the first term in Equation (HIT) is negli¬ 
gible, and so the condition for E{Rg) > 0 can be sum¬ 
marized as 


gdm Mg 

Eg \M\ 


> - (3 - 2jEg) 

7] ^ ^ 


(43) 


or (using Equations (|37|) and ([38|) ) 


tadv^G 

^heat,G 


> - (3 - 2^Eg) ■ 
V 


(44) 


For the typical values of 77 = 0.6, ^Eg = 1-4 and zero dis¬ 
sociation, the expression on the right hand side is about 
0.3. A comparison of this result to actual simulations 
is presented in Figure 1131 Indeed, runaway expansion 
is seen to commence in conjunction with the ratio of 
timescales in the gain region exceeding this threshold. 

We note that dissociation losses can be included in the 
estimate for the required value of tadvG/^heat,G by cor¬ 
recting the non-adiabatic integrand to Q{r) — qd- While 
this approximation is somewhat crude (since the quasi- 
stationary approximation is less appropriate for large dis¬ 
sociation losses), it does provide some quantitative esti¬ 
mate of the effect. The actual derivation then yields a 
correction of +{qdMG)/[E gvi) to the right hand of Equa¬ 
tion (1441) . demonstrating that a larger dissociation energy 
requires a larger timescale ratio, since heating must com¬ 
pensate for the lost energy. With the aid of the quasi- 
stationary approximation we predict that in the case of 
full dissociation the necessary value of tadvG/theat,G is 
about 1.4. Again, this value appears to be consistent 
with the simulations, as demonstrated in Figure [131 

Finally, we point out that the value of tadv,G/theat,G 
increases to the from left to right (increasing radii) in the 
Rs — Vs plane, so a ratio of unity should in general not 
be far from the critical 05 = 0 curve. However, the two 
are by no means identical criteria, and we find no distinct 
relation between tadv,G/theat,G and any specific feature 
in the phase space. Moreover, the value of unity does 
not appear to hold any particular significance: it occurs 
well into runaway in the absence of dissociation, while it 
is an insufficient condition (occurs during the oscillations 
prior to runaway) when full dissociation is assumed. 

We summarize our analysis of the ratio between advec- 
tion and heating time scales with two claims: 


• Strictly speaking, tadv,G/theat,G > 1 is not a special 
condition for an explosion. Rather, it is a value 
that apparently lies in the path of a shock that 
is on its way to a runaway expansion, and that it 
occurs close to (before or after) the transition to a 
positive shock acceleratioifl 


where Q{r) was defined in Equation (QH). Since g > 0 in 
the gain region, Q(r) monotonically decreases with the 
radius, and Q{Rs) = 0 . The total non-adiabatic integral 
over the entire gain region can therefore be approximated 
by: 

(•Rs 


( 3 (r)dm = r]Q{RG)MG 


(42) 


Rg 


• We also do not identify the gain region as being 
in any way unique in terms of a necessary condi¬ 
tion for runaway expansion. The curves in Figure 

^ A similar argum ent was raised independently by 
IMurphv He Dolencel II2015I ). just as this manuscript was being 
completed. 
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Figure 11. Phase-space diagram (see Figure !^ with constant E = 0 lines. See text for details. 



Figure 12. The total effective energy E(r) for |Mo| = O.SMq s“^, 
L 52 = 8 and = 0 for different Rs{t = 0). The locations of Rq 
and Re are marked by circles and squares, respectively. See text 
for definitions of these quantities. 
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Figure 13. The ratio of the advection timescale to heating 
timescale as a function of time in simulations with |Mo| = 
0.8 Mq s“^, and different neutrino luminosities L, in units of 
lO^^erg s~^, Rs,o = 100 km. Solid lines with = 0 and dashed 
lines with q^ = 8.5 x lO^^erg g~^. The ratios should be com¬ 
pared to the threshold of tadv cl'^heat G — 0.3, predicted by the 
quasi-stationary approximation for = 0 (black solid line) and 
tadv,Gliheat,G = 1-4, predicted for for qd = 0. 


imply that M{Re) > 0, or some other alterna¬ 
tive curve, could be used just as well a necessary, 
but not sufficient, condition. Clearly, runaway ex¬ 
pansion requires signihcant heating, but evolution 
towards a runaway expansion is a quality of the 
entire accretion layer, rather than just of the gain 
region. 

7. SUMMARY AND DISCUSSION 

In this work we investigated the driving of a stalled 
accretion shock into runaway expansion in the context 
of the neutrino heating mechanism in core-collapse su¬ 
pernovae. We used spherically symmetric simulations in 
conjunction with analytic derivations, and examined the 
evolutionary paths the system may take towards runaway 
expansion, if it occurs. Our numerical results included 
a parameter survey of the neutrino luminosity and the 
incoming accretion rate, as well as of the shock parame¬ 
ters, namely, its instantaneous radius and velocity. Our 
main conclusions are as follows. 

1. The instantaneous shock acceleration, as, depends 
on both the shock radius, Rs, and velocity, Vs- 
We examine the shock acceleration in the Rs — Vs 
plane, and show that it is generally divided into re¬ 
gions of positive and negative shock acceleration 
(see Figures |4| and |5|). We expanded the com¬ 
monly used stationary approximation, by includ¬ 
ing a non-zero shock velocity and solving the re¬ 
sulting profile in the accretion layer with modified 
boundary conditions I TS.II) . This quasi-stationary 
approximation developed here allows to map this 
plane in advance of a time-dependent simulation, 
and nicely predicts the evolution of a model. The 
quasi-stationary approximation is especially accu¬ 
rate when the shock dynamical time, ts = Rs/^S, 
and oscillation period, tosc, are much longer than 
the advection time through the accretion layer, 
which is the case if the immediate dissociation en¬ 
ergy losses across the shock are not too large. We 
also use the derivation to estimate the oscillation 
period, and show why it is typically 25-50 ms, very 
weakly dependent on the mass accretion rate (see 
















































Phase Space Analysis of Core-Collapse Supernova 


17 


Appendix [C|). This result is consistent with simu¬ 
lations. 

2. Generally, the Rs — Vs phase space shows positive 
shock acceleration at small and large shock radii, 
with an intermediate trough of negative accelera¬ 
tion. The trough is bounded by two as = 0 curves. 
The one corresponding to smaller radii is an ’’os¬ 
cillatory” as = 0 curve: if the shock moves from it 
inward, positive acceleration due to the high pres¬ 
sure close to the PNS can cause the shock to bounce 
back and oscillate. Such oscillations either damp 
out or grow, depending on the parameters of the 
flow. In contrast, the as = 0 curve at large radii is 
” critical”, since crossing it with a positive velocity, 
leads to a region of positive acceleration, so run¬ 
away expansion is guaranteed (i.e., profiles which 
correspond to this curve are inherently unstable). 
The magnitude of the neutrino luminosity comes 
into play through quantitative aspects of this phase 
space structure. For low luminosities, the trough 
of negative acceleration is both deep and wide, and 
when the initial shock radius and velocity corre¬ 
spond to this trough, this shock is generally stable 
and oscillations are damped. When the neutrino 
luminosity is high, the trough of negative accel¬ 
eration becomes shallow and narrow, and larger 
regions of it correspond to profiles which are un¬ 
stable, so the shock eventually evolves beyond the 
critical as = 0 curve and towards runaway expan¬ 
sion. 

3. For arbitrary combinations of Rs and Vs serving as 
initial conditions, the evolution of a model can be 
tracked in the phase space estimated by the quasi¬ 
stationary approximation. Our analysis demon¬ 
strates that the shock velocity must be accounted 
for when assessing whether the evolution will end 
in runaway expansion or not. 

4. In the special case of initial stationary profiles 
(Vs = 0, as = 0), we find that the partial deriva¬ 
tives of the shock acceleration offer a satisfactory 
indication regarding the outcome. If das/dRs > 0 
and das/dVs < 0 the initial profile lies on the 
critical as = 0 curve, and the profile will evolve 
to runaway expansion in a non-oscillatory fashion. 
In contrast, if das/dRs < 0 the initial profile lies 
on the oscillatory as = 0 curve and its evolution 
depends on the sign of das/dVs- If this deriva¬ 
tive is also negative, the system is stable and any 
small initial perturbation will damp out and set¬ 
tle to the stationary profile. On the other hand, 
if das/dVs > 0 , oscillations will tend to grow un¬ 
stably (be ”anti-damped”), leading to runaway ex¬ 
pansion through an oscillatory path. This analy¬ 
sis bares out very well in the simulations, in the 
sense that we can predict the critical luminosity 
as a function of mass accretion rate, Lc(Mq), for 
both oscillatory and non-oscillatory explosions; see 
Figure [T] Quantitatively, das / dVs = 0 is a very 
good indicator for finding the critical luminosity 
for oscillatory explosions, and das/dRs = 0 does 
well in predicting the critical luminosity for non- 
oscillatory explosions. We note that we slightly 


underestimate Lc(Mo) for the oscillatory mode, be¬ 
cause there is a small range of luminosities which 
generate stable oscillations. 

5. Both modes of runaway expansion are available for 
low mass accretion rates. In this case, oscillatory 
runaway with das/dVs > 0 is met at lower lumi¬ 
nosities than das/dRs > 0 , so the actual critical 
luminosity is due to anti-damped oscillations. We 
find that the luminosity at which das/dVs = 0 , 
Lc,osc, is strongly dependent on the mass accre¬ 
tion rate, and there exist some finite rate Me for 
which Tc,osc diverges. At larger accretion rates, 
\Mo\ > Me-, only non-oscillatory runaway expan¬ 
sion is possible, with a critical luminosity Lc^dir-, 
which is weakly dependent on Mq. The value of 
Me may depend on the specifics of the model (PNS 
properties, initial shock radius, dissociation losses, 
etc.), but for reasonable parameters appears to be 
about 1 Mqs~^. 

6 . We also applied the quasi-stationary approxima¬ 
tion to examine the commonly-used ratio between 
the advection time scale and the heating time scale 
in the gain region of the accretion layer (see 

We find that for the gain region to have a posi¬ 
tive effective energy, this ratio must be at least a 
few tenths, with some dependence on the specific 
loss of energy due to dissociation across the shock 
(Figure [131). This result implies that the ratio of 
time scales, while not a fundamental condition for 
an explosion, can serve as an indication if a given 
system will evolve to an explosion, as suggested in 
previous works. However, we note that there are 
several qualifications to this observation. First, a 
positive energy in the gain region appears to be a 
threshold that must be crossed on the path to run¬ 
away expansion, yet it is not an actual condition 
for an explosion, since it does not exactly coincide 
with the region of positive shock acceleration (Fig¬ 
ure El). Second, we find that in many cases the ad¬ 
vection time becomes longer than the heating time 
only after the transition to a runaway expansion 
has occurred, implying that an equality between 
the two time scales can be a consequence of a suc¬ 
cessful explosion, rather than a physical condition 
for initiating one. Third, we do not identify the 
gain region as being unique in some physical sense; 
in principle, a alternative criteria can be formulated 
by applying it to a different region in the accretion 
layer, such as that which holds the largest total en¬ 
ergy. We therefore suggest that the properties of 
the gain region should not be used as a singular 
measure of the likelihood of an explosion. 

Since our simulations are inherently limited to spher¬ 
ical symmetry (with additional assumptions and simpli¬ 
fications), quantitative differences are to be expected in 
the actual problem of the neutrino driven mechanism. 
Nonetheless, we do believe that the qualitative nature 
of our results is applicable to the full, multi-dimensional 
problem. In particular, we emphasize the dependence of 
the shock acceleration on its radius and velocity; this fact 
certainly facilitates dynamical spherical simulations to 
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explode with a lower neutrino luminosity than predicted 
by the stationary approximation. We speculate that this 
behavior does carry over to the three dimensional reality, 
where the shock can sample a wider range of combina¬ 
tions of Rs and Vs than in the one dimensional case. 
This should result from initial aspherical co nditions in 
the silicon-oxygen layer (iCouch fc OttI [2M^ . and also 
from the onset of turbulence (of course, turbulent pres¬ 
sure should_alsocOTitr^ute_to_^£^tatiM an explosion, 
fsee lCouch fc OttI (j201,^ : lFernanded ()2015h 'll. This com¬ 
bination may be responsible - at least in part - for the 
observation that the critical luminosity in two and three 
dimensional simulations is lower than in the one dimen¬ 
sional case. 

A natural expansion of this work would be to examine 
two- and three dimensional simulations, and to quantita¬ 


tively compare the actual acceleration of the shock with 
that predicted from a spherical model of identical pa¬ 
rameters. Such a comparison will clarify the dependence 
of the shock characteristics upon the physical processes 
inside the accretion layer, such as turbulence and con¬ 
vection. 
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APPENDIX 

A. THE DERIVATION OF (iPl/dt^)Qs 

In this appendix we derive the various terms which appear in the quasi-stationary approximation for <Pl/dt^ (equa¬ 
tion [531) ■ Recall that the goal of the approximation is to include the effect of a finite shock velocity on the profile of 
the accretion layer. To do so, we assume that the profile adjusts instantaneously to changes in the shock radius and 
velocity. This requires some care, as we show below, because the advection time is not negligible. 


A.l. Virial theorem for spherical accretion 

In general form, the first and second time derivatives of the moment of inertia in the accretion layer are: 


<U 

dt 


Sirr^pudi — 


r'^Mir) — d-KT^pV 


1 S 


J P 


where we used the continuity equation (l2|), and 


d^I 
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o s dipu) , 
87 rr ^——-dr — 
dt 


IRP 


— — (dTrr^pU) 


(Al) 


(A2) 
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where M{r) = Arir'^pu is the local mass accretion rate, M = dM/dt is its time derivative, and V is the velocity of 
the flow boundaries. As was defined in (JH the notation [• • • ]p stands for the difference bet ween the expression in the 
immediate downstream of the shock and just outside the PNS. The last term in Equation (IA2I) is a Lagrangian time 
derivative of the boundaries. 

Using the momentum conservation equation ([3]) we express dfl/dt^ in a ’’virial theorem” fashion with boundary 
conditions: 

1 d^I 

+ ^^PNS + AEs 
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and K and are the total kinetic and gravitational energies, defined as 
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The physical interpretation of E and AEpjys + AEs is the energetic contributions of the accretion layer and of the 
boundaries (at the PNS and at the shock, respectively). We note that the boundary terms are calculated within the 
accretio n layer, just inside the shock at Rs and just above the PNS at Rp. The time derivative d/dt in Equations 
(IA5IIA6I) is a Lagrangian derivative of values in the post-shocked region, which is identically zero for a stationary PNS 
surface. _ 

The dl^/dt^ expression in Eqs. (IA3MA8I) is precise, but does not give accurate results in the context of a quasi-steady 
model. This occurs because the P/p term in E is everywhere distorted by the motion of the shock, which is nonphysical 
when the advection time is not negligible with respect to the characteristic time for changes in the shock velocity, e.g., 
the oscillation time. Therefore, we next seek an alternative formulation that is more susceptible for the quasi-steady 
approximation. 


A.2. The Effective Energy Term 

We now turn to estimate the effective energy term E in Equation (1231) . using the quasi-stationary approximation 
by decomposing it into several terms. Integrating the energy-conservation equation (l4|) over volume in a stationary 
profile {d/dt = 0) profile yields: 

M[bs-b{r)] = Q{r), (A9) 

where b{r) is the Bernoulli function (specific energy) at r, 


j.i' ^ ^ 2 I , P GMp 

b(r) = -u + e-\ - - 

2 p r 


(AlO) 


bs = b{Rs), and Q{r), also defined in the main text (Equation (|29)l '). is the total non-adiabatic energy contribution 
rate between r to Rs, 

pRs pRs 

Q{r) = / dTrpr^qdr = / gdm. (All) 

Jr Jr 

Using equation ([9l) for energy conservation across the shock yields 


bs = bQ- qd + Vs{ui - uo) = -Qd - 




(A12) 


where the Bernoulli function just outside the shock, 69 = 0, vanishes in the pre-shocked material (Equation!^. The 
Bernoulli function inside the post-shocked flow is then obtained by 


1 ! \ 1 2 P GMp 1 ^/ \ 7 

b{r) = -u^ + e+ -= ^Q{r) + bs , 

2 p r \M\ 


and the ratio P/p in the shock can be rewritten as 


P 

P 


7-1 

7 


\M\ 


Q{r) + bs 


GMp 

r 



(A13) 


(AM) 


where 7 is the local adiabatic index of the flow. The total non-adiabatic energy gained by the material in the accretion 
layer is then 

U = [ 77 ^‘5(0dm . (A15) 

Jrp \M\ 

Substituting these derived quantities into the integrals included in the effective energy E from Equation (I23p . we 
arrive at a form for the effective total energy in the quasi-stationary approximation, Eq 5 


EQs=K + U + n + Bs (A16) 

The terms in Equation (IA16I) are effective manipulations of the original functions, including coefficients which depend 
on the local adiabatic index in the accretion layer. We define K, U, (l and Bs as the effective kinetic energy, effective 
gained energy, effective gravitational potential and effective Bernoulli function, respectively: 


K = 


I R„ 


3 - 7 ^ 1 2 7 
—u dm 


n = 


3-27 


GMr, 


dm 


(A17) 

(A18) 




(A19) 
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and 




bs{m) dm . 


(A20) 


The last expression requires some discussion. The quantity bs in the integral should reflect the Bernoulli function 
bs(jn) of the mass element m when it crossed the shock, and so is not a constant across the accretion layer. For example, 
when dissociation is neglected, bs = Vs{us — uq), which vanishes in the stationary case, and fluctuates between positive 
and negative values as the shock oscillates. In general, this causes the integral Bs to be largely negligible with respect 
to the other terms of the effective energy, which we confirm explicitly in the simulations. Hence, in the quasi-stationary 
approximation: 

Eqs ^K + i/ + n. (A21) 


A.3. The Boundary Terms 

The boundary terms in d^I/dt^ includes distinct contribution from the PNS and from the shock. A key feature 
of our model is that the boundary contributions AEpvrs and AEs (Eqs. (IA5HA6I1 ') are invariant under the Rankine- 
Hugoniot relations Eqs. (171191) . Therefore, we can determine the boundary terms using the known conditions outside 
the postshocked region between the PNS and the shock. 

At the PNS radius the inner boundary can be estimated with the quantities just below Rp, which include the 
predetermined pressure, Pp. The last term in equation (IA5I1 is null since Rp is fixed {Vp = 0), and so 

AEpjvs = 47ri?p (Pp + Pf Up) + —Ppfifp . (A22) 

The result in Equation (IA22I) is general, but typically thermal pressure dominates over ram pressure at the PNS and 
we can approximate Pp + ppu^p ~ Pp. In fact, in constructing the phase space, we found that accuracy can be 
increased further by imposing the fixed total pressure at the Pp as the sum of the actual thermal pressure at the PNS 
and the initial value of PpUp, which is a correction of a few percents. Thus the approximation is only that PpUp is not 
updated as the flow evolves, and this error is limited to less than one percent. In the quasi-stationary approximation 
M is uniform in space (although allowed to vary in time). Hence, we set Mp M and Mp cs M and use equation 

m- 

At the shock, the quantity AE5 is also invariant under the Rankine-Hugoniot relations can be calculated from the 
upstream values (index 0) instead of the shocked values (index I): 

AEs = (Pi + psuj) + ip|M(P5) - 4 (27rP|p«ys) 

1 .. d 

= 47rP| (Pq -|- PoUq) 4" ~ ^ (27rP|poRs) 

Einally, for accretion through the shock arriving as pressure-less free-fall at a constant accretion rate, Mq (so Mq = 0), 
and Ws is equal to: 

AE. = a(GMpy^im («!'“)) (A24) 

Por further application in the quasi-stationary model we denote the work done by the PNS and the energy advected 
across the boundaries by Wpns and Wp, respectively, as 

WpNS - ‘ii^RpPp , (A25) 


Wb ^ -a(GMp)i/2|Mo|Py^ -h 


|Mo| d 

7a(GMp)i/2 dt 



(A26) 


The equation is approximate due to assuming a uniform accretion rate and neglecting the ram pressure. Note that 
under the above assumptions we maintain the sum Wpns + Wp — AEpjvs -|- AE5. Note that 5 = (/3 — 1) {Rp/RsY 
being at the order of unity. 

It is noteworthy that the shock related terms do include an explicit dependence on the shock radius. The first term 
in Equation (|A26p is negative and restr ains the shock, corresponding to the impulse of the accreting matter on the 
shock, while the term in Equation (|A25ll is positive (assisting shock expansion) and accounts for the kinetic energy of 
the shocked matter and the work it does on the shock. 

Summarizing the results we conclude that the second time derivative of the spherical moment of inertia in the 
quasi-stationary approximation presented in the main text is given by: 


1 d^J 

2 dt^ 


k + u + n + 47rP|,(Pp -h ppu],) - a(GMp)i/2|Mo|Py^ + 


|Mo| d 

7a{GMpy/^ dt 



(A27) 
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B. ON THE IMPORTANCE OF THE DIEFERENT TERMS IN EQUATION JSS]) 

Our final estimate for the shock acceleration in Equation (IM)l includes the effective energy of the accretion layer, 
]Eq 5 , the work done by the PNS, Wpjvs, and two terms which depend on the properties of the shock. Here we show 
that WpNS^ which is invariably positive, generally dominates over the shock terms. Correspondingly, a positive Eqs 
can serve as a sufficient (albeit not necessary) condition for a positive shock acceleration. 

First, recall that Pp is calculated through a stable, stationary accretion flow which corresponds to a point on the 
V 5 = 0 axis and is bound such that Eqs < 0. Requiring cPl/dt^ = 0 along with Egs < 0 in equation (|23l) implies that 
in the stationary solution, the pressure term at the PNS must dominate over the mass influx at the shock: 

A-kR%Pp > a(GMp)i/2|Mo|Ry^ = \MoRsuo\ (Bl) 

This inequality holds as long as the profile is still in the oscillatory regime. It becomes invalid once the flow has reached 
runaway expansion (large shock radii), but by then the explosion is assumed to be well under way with a large shock 
velocity, and our quasi-stationary approximation breaks down in any case. 

The second term originating from the shock in Equation (1551) is {—Vgd(/dRs), with ( defined in Equation (1511) . 
Taking the partial derivative with respect to the shock radius yields 


..2 iP-^)\Mo\ 

^ dRs ^ 2a{GMpy/^ 


[Rg - Rp) 



+ 2RsRs^^^ 


f5 

ul 2 ^2 




(B2) 


The equality is approximate because we nefglected the weak dependence of j3 on the shock radius. Now, for the 
small shock velocities we consider in the oscillatory phase, Vs/uq << 0 . 1 , the coefficient of \MoRsuo\ in the expression 
in Equation (IB2I) is much smaller than unity (/3 ~ 6 — 10 and (5/2 — i?p/2i?|) < 5/2). Correspondingly, Vgd(/dRs 
is subdominant to the first shock related term (Equation (IBII) '). and does not change the hierarchy in which the PNS 
term dominates. 


C. SHOCK OSCILLATIONS TIMESCALES 


In this appendix we show that the virial theorem can be used to explain the typical oscillation period of tens 
of milliseconds we find in the simulations. Coming back to Equations (1T51) . (IA3I) and (IA22I - [AM1) we can write a 
oscillator-like equation for the shock radius when assuming a uniform accretion rate in the postshocked region, 

^ {2TT^PoRtVs) ~ E + inR^pPp - \MoRsuo\ (Cl) 

where (j) = {/3 — 1)[1 — (Rp/Rs^] and E is define d in equation IA4I (i.e., the quasi-stationary approximation for the 
energy is not necessary in this context). Equation (IC1() was found to be in good agreement with the shock motion for 
both zero and full dissociation parameters. For small oscillations around a stationary solution for which Vg = 0 at 
Rs, 0 j the right hand side must be zero at Rs,o- 

We now assume that during the oscillations, the change in the effective energy is roughly proportional to the change 
of the inertia crossing the shock: AE = p\MoRsuo\, with 0 < ^ < 1. This is a lowest-order approximation, in which 
the change in E depends only on the shock radius (the dependence on the shock velocity is neglected). We confirmed 
this assumption quantitatively in the simulations. It conveys the fact that in small oscillations the accretion layer 
adjusts to include the material that is either added or lost as the shock moves. Correspondingly, (recalling that Wpns 
is constant in the quasi-stationary approximation): 


d r (j) 

dt [ 2 a(GMp)i /2 


I Mo 1 4 /Vs 


-(1 - p)a{GMpy/^\My (i?4 - i?4) 


(C2) 


Considering that the shock radius is significantly larger than the PNS radius, 1 — (Rp/RsY ~ 1 and combined with 
the fact that the compression ratio depends weakly on the shock radius, we can treat (j) as roughly constant during 
the oscillations. Equation (IC2I) can now be rewritten as 


dV4 7(1 - M)aVMp / ^ 1/2 „i /2 

dP ~ 6 VV V,o 


(C3) 


It is noteworthy that the mass accretion rate Mq has canceled out of the equation. For small oscillations, this is an 
harmonic oscillator equation with a time period of: 


To 


|mo| y 1 - m 


(C4) 


For a compression ratio in the range /3 = 6 — 10 and Rp/Rs <0.5, the value of is confined to values of 2-3. 

Unless p. is very close to unity, we conclude that the oscillation period should be (10 — 20) i?s, o/|mo|, or (25 -50) 

milliseconds. This is indeed in good agreement with the results of the simulations. We also recover the general trend 
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that the period should be weakly dependent on the mass accretion rate and on the neutrino luminosity, as their effect 
is limited to the finer details of \/^/(l — ^). 

Finally, we note that in the approximation above, the oscillations are unconditionally stable. This is due to the 
neglect of a 9(AE)/9Vs term in the derivation. This partial derivative is directly related to the das/dVs derivative 
discussed in the main text, which determines whether the shock oscillations around the stationary solutions will damp 
or grow to a runaway expansion. 
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